

/*
 * Class:        NormalInverseGaussianProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;
import umontreal.iro.lecuyer.hups.*;



/**
 * This class represents a normal inverse gaussian process (NIG).
 * It obeys the stochastic differential equation
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="eq:nig"></A>
 * <I>dX</I>(<I>t</I>) = <I>&#956;dt</I> + <I>dB</I>(<I>h</I>(<I>t</I>)),
 * </DIV><P></P>
 * where 
 * <SPAN CLASS="MATH">{<I>B</I>(<I>t</I>),&nbsp;<I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN> is a {@link BrownianMotion} 
 * with drift <SPAN CLASS="MATH"><I>&#946;</I></SPAN> and variance 1,
 * and <SPAN CLASS="MATH"><I>h</I>(<I>t</I>)</SPAN> is an {@link InverseGaussianProcess} 
 * <SPAN CLASS="MATH"><I>IG</I>(<I>&#957;</I>/<I>&#947;</I>, <I>&#957;</I><SUP>2</SUP>)</SPAN>, with 
 * 
 * <SPAN CLASS="MATH"><I>&#957;</I> = <I>&#948;dt</I></SPAN> and 
 * <SPAN CLASS="MATH"><I>&#947;</I> = (&alpha;^2 - &beta;^2)<SUP>1/2</SUP></SPAN>.
 * 
 * <P>
 * In this class, the process is generated using the sequential
 * technique:  <SPAN CLASS="MATH"><I>X</I>(0) = <I>x</I><SUB>0</SUB></SPAN> and 
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay">
 * <I>X</I>(<I>t</I><SUB>j</SUB>) - <I>X</I>(<I>t</I><SUB>j-1</SUB>) = <I>&#956;dt</I> + <I>&#946;Y</I><SUB>j</SUB> + (Y_j)<SUP>1/2</SUP><I>Z</I><SUB>j</SUB>,
 * </DIV><P></P>
 * where 
 * <SPAN CLASS="MATH"><I>Z</I><SUB>j</SUB>&#8764;<I>N</I>(0, 1)</SPAN>, and 
 * <SPAN CLASS="MATH"><I>Y</I><SUB>j</SUB>&#8764;<I>IG</I>(<I>&#957;</I>/<I>&#947;</I>, <I>&#957;</I><SUP>2</SUP>)</SPAN>
 * with 
 * <SPAN CLASS="MATH"><I>&#957;</I> = <I>&#948;</I>(<I>t</I><SUB>j</SUB> - <I>t</I><SUB>j-1</SUB>)</SPAN>.   
 * 
 * <P>
 * There is one {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} 
 * used to generate the <SPAN CLASS="MATH"><I>Z</I><SUB>j</SUB></SPAN>'s and
 * there are one or two streams used to generate the underlying
 * {@link InverseGaussianProcess}, depending on which IG subclass
 * is used.
 * 
 * <P>
 * In finance, a NIG process usually means that
 * the log-return is given by a NIG process;
 * {@link GeometricNormalInverseGaussianProcess}
 * should be used in that case.
 * 
 */
public class NormalInverseGaussianProcess extends StochasticProcess  {
    protected RandomStream streamIG1;   // U[0,1) gen used in the inverse gaussian
    protected RandomStream streamIG2;   // U[0,1) gen used (maybe) in the inverse gaussian
    protected RandomStream streamBrownian;   // U[0,1) gen for the normal "Brownian"

    protected InverseGaussianProcess igProcess; 
    protected NormalGen normalGen;
    // Randomized time generated by the InverseGaussianProcess.
    protected double[] stochTime;
    double[] dt;
    double[] mudt;

    // Parameters of the normal inverse gaussian
    protected double mu;
    protected double delta;
    protected double alpha;
    protected double beta;
    protected double gamma;





   /**
    * Given an {@link InverseGaussianProcess} <TT>igP</TT>, constructs a
    * new <TT>NormalInverseGaussianProcess</TT>.  The parameters and observation
    * times of the IG process will be overriden by the parameters
    * of the NIG process.  If there are two 
    * {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}'s in the 
    * {@link InverseGaussianProcess}, this constructor assumes that
    * both streams have been set to the same stream.
    * 
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        InverseGaussianProcess igP)  {
        setParams (x0, alpha, beta, mu, delta);
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)
        igProcess = igP;  // its params will be set in init().
        this.streamIG1 = igProcess.getStream();
        this.streamIG2 = streamIG1;
    }


   /**
    * Constructs a new <TT>NormalInverseGaussianProcess</TT>.
    * The string argument corresponds to the type of underlying
    * {@link InverseGaussianProcess}.  The choices are SEQUENTIAL_SLOW,
    * SEQUENTIAL_MSH, BRIDGE and PCA, which correspond
    * respectively to {@link InverseGaussianProcess}, {@link InverseGaussianProcessMSH},
    * {@link InverseGaussianProcessBridge} and {@link InverseGaussianProcessPCA}.
    * The third {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}, <TT>streamIG2</TT>,
    * will not be used at all if the SEQUENTIAL_SLOW or PCA methods are chosen.
    * 
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        RandomStream streamIG1,
                                        RandomStream streamIG2,
                                        String igType)  {
        setParams (x0, alpha, beta, mu, delta);
        this.streamIG1 = streamIG1;
        this.streamIG2 = streamIG2;
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)

        // The initial time of igProcess here, 0., is arbitrary: 
        // init() sets igProcess.x0 = t[0].
        if(igType.compareTo("SEQUENTIAL_SLOW") == 0)
            // the initial value, 0.0, of the IG will be overriden by
            // the initial time, when it will be set.
            igProcess = new InverseGaussianProcess(0.0, delta, gamma, 
                                streamIG1);
        else if(igType.compareTo("SEQUENTIAL_MSH") == 0)
            igProcess = new InverseGaussianProcessMSH (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("BRIDGE") == 0)
            igProcess = new InverseGaussianProcessBridge (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("PCA") == 0)
            igProcess = new InverseGaussianProcessPCA (0.0, delta, gamma, 
                                streamIG1);
        else throw new IllegalArgumentException("Unrecognized igType");
    }


   /**
    * Same as above, but all {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}'s 
    * are set to the same stream, <TT>streamAll</TT>.
    * 
    */
   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamAll,
                                        String igType)  {
        this(x0, alpha, beta, mu, delta, streamAll, streamAll, streamAll, igType);
    }


   /**
    * Generates the path.  
    * This method samples each stream alternatively,
    * which is useful for quasi-Monte Carlo, where
    * all streams are in fact the same iterator on 
    *  a {@link umontreal.iro.lecuyer.hups.PointSet PointSet}.
    * 
    */
   public double[] generatePath()  {
        if(igProcess.getNumberOfRandomStreams() != 1)
            return generatePathTwoIGStreams();

        double x = x0;
        double[] randomNormal = new double[d];
        double[] randomIG1 = new double[d];
        for (int j = 0; j < d; j++)
        {
            randomIG1[j]    = streamIG1.nextDouble();
            randomNormal[j] = streamBrownian.nextDouble();
        }

        stochTime = igProcess.generatePath(randomIG1);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(randomNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }


    protected double[] generatePathTwoIGStreams() {
        double x = x0;
        double[] uniformNormal = new double[d];
        double[] uniformIG1 = new double[d];
        double[] uniformIG2 = new double[d];

        for (int j = 0; j < d; j++)
        {
            uniformIG1[j]    = streamIG1.nextDouble();
            uniformNormal[j] = streamBrownian.nextDouble();
            uniformIG2[j]    = streamIG2.nextDouble();
        }

        stochTime = igProcess.generatePath(uniformIG1, uniformIG2);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(uniformNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }

   /**
    * Returns the value of the process for the next time step.
    * If the underlying {@link InverseGaussianProcess}
    * is of type {@link InverseGaussianProcessPCA}, this method cannot be
    * used.  It will work with {@link InverseGaussianProcessBridge}, but
    * the return order of the observations is the bridge order.
    * 
    */
   public double nextObservation()  {
       double igNext = igProcess.nextObservation();
       observationIndex = igProcess.getCurrentObservationIndex();
       stochTime[observationIndex]  = igNext;
       double dY = igNext - stochTime[0];
       double dT = t[observationIndex] - t[0];
       path[observationIndex] = x0  +  mu * dT  +  beta * dY  +
                                Math.sqrt(dY) * normalGen.nextDouble();
       return path[observationIndex];
    }


    protected void init()
    {
        super.init();
        igProcess.setParams(delta, gamma);
        if(observationTimesSet){
            stochTime = new double[d+1];
            stochTime[0] = t[0];
            dt = new double[d];
            mudt = new double[d];
            for(int i = 0; i < d; i++){
               dt[i]   = t[i+1] - t[i];
               mudt[i] = dt[i] * mu;
            }
        }
    }

   /**
    * Sets the observation times on the NIG process as usual,
    * but also sets the observation times of the underlying {@link InverseGaussianProcess}.
    * It furthermore sets the starting <SPAN  CLASS="textit">value</SPAN> of the {@link InverseGaussianProcess}
    * to <TT>t[0]</TT>.
    * 
    */
   public void setObservationTimes(double t[], int d)  {
        super.setObservationTimes(t, d);
        igProcess.setObservationTimes(t, d);
        igProcess.x0 = t[0];
    }


   /**
    * Sets the parameters. Also, computes 
    * 
    * <SPAN CLASS="MATH"><I>&#947;</I> = (&alpha;^2-&beta;^2)<SUP>1/2</SUP></SPAN>.
    * 
    */
   public void setParams (double x0, double alpha, double beta,
                          double mu, double delta)  {
        // Initializes the parameters of the process
        if (delta <= 0.0)
            throw new IllegalArgumentException ("delta <= 0");
        if (alpha <= 0.0)
            throw new IllegalArgumentException ("alpha <= 0");
        if (Math.abs(beta) >= alpha)
            throw new IllegalArgumentException ("|beta| >= alpha");

        this.x0    = x0;
        this.mu    = mu;
        this.delta = delta;
        this.beta  = beta;
        this.alpha = alpha;

        gamma = Math.sqrt(alpha*alpha - beta*beta);
        if (observationTimesSet) init();
    }


   /**
    * Returns alpha.
    * 
    */
   public double getAlpha()  {
      return alpha;
   }


   /**
    * Returns beta.
    * 
    */
   public double getBeta()  {
      return beta;
   }


   /**
    * Returns mu.
    * 
    */
   public double getMu()  {
      return mu;
   }


   /**
    * Returns delta.
    * 
    */
   public double getDelta()  {
      return delta;
   }


   /**
    * Returns gamma.
    * 
    */
   public double getGamma()  {
      return gamma;
   }


   /**
    * Returns the analytic average, which
    * is 
    * <SPAN CLASS="MATH"><I>&#956;t</I> + <I>&#948;t&#946;</I>/<I>&#947;</I></SPAN>.
    * 
    */
   public double getAnalyticAverage (double time)  {
        return mu*time+delta*time*beta/gamma;
    }


   /**
    * Returns the analytic variance, which is
    * 
    * <SPAN CLASS="MATH"><I>&#948;t&#945;</I><SUP>2</SUP>/<I>&#947;</I><SUP>3</SUP></SPAN>.
    * 
    */
   public double getAnalyticVariance (double time)  {
        return delta*time*alpha*alpha/gamma/gamma/gamma;
    }


   /**
    * Only returns the stream if all streams are equal,
    * including the stream(s) in the underlying {@link InverseGaussianProcess}.
    * 
    */
   public RandomStream getStream()  {
        if( (streamIG1 != streamIG2)                ||
            (streamIG1 != streamBrownian)         ||
            (streamIG1 != normalGen.getStream())  ||
            (streamIG1 != igProcess.getStream())  )
            throw new UnsupportedOperationException
            ("Two different streams or more are present");
        return streamIG1;
    }


   /**
    * Sets all internal streams to <TT>stream</TT>,
    * including the stream(s) of the underlying {@link InverseGaussianProcess}.
    * 
    */
   public void setStream (RandomStream stream)  {
        streamIG1 = streamIG2 = streamBrownian = stream;
        normalGen.setStream(stream);
        igProcess.setStream(stream);
    }


} 

